METHOD FOR PACKING AND FILTERING GEOPHYSICAL EVENTS TO 
RETRIEVE THEREFROM DATA ON THE TYPE OF SUBSOIL 



BACKGROUND OF THE INVENTION 
Field of the Invention 

[0001] The present invention relates to a method for packing and filtering 

geophysical data in order to retrieve therefrom information on the nature of the subsoil. 

Description of the Prior Art 
[0002] Geophysical measurements are conventionally used to provide 

complementary information, in relation to direct observations in wells, on the geologic 
properties variations of underground formations. Seismic interpretation is for example 
based on the analysis of attributes in order to retrieve these geologic properties from the 
seismic information, in a zone of interest of an underground formation (reservoir zone for 
example). Another example is log interpretation where the various measurements are 
simultaneously analyzed to deduce therefrom the petrophysical characteristics of the 
underground formation studied. 

[0003] A standard approach for seismic or log attributes analysis involves the use 

of statistical pattern recognition and estimation methods in order to calibrate the pertinent 
attributes with the available information obtained in various wellbores. 
[0004] For example, French Patent 2,768,818 and corresponding US Patent 

6,051,651 filed by the assignee describes a method allowing detection of the potential 
classes in a population of events related to the physical properties of a complex medium 
such as the subsoil, located from data obtained by subsoil exploration or by in-situ 



measurements (events of geologic, geophysical nature, events linked with fluid production 
data, etc.), these events being located by points in a multi-variable space and defined by an 
attributes vector. 

[0005] Another example is given by French Patent Application 01/05,675 filed by 

the assignee, which describes a method for facilitating monitoring over time of the 
evolution of an underground zone by compared analysis of a certain number n of seismic 
records successively obtained after n successive 3D seismic surveys (4D seismic method), 
wherein a pattern recognition technique applied to the whole of the seismic events of 
several surveys, considered and analyzed simultaneously, is used. 

[0006] EP Patent A-671,017 and corresponding US Patent 5,638,269 filed by the 

assignee can also be mentioned, which describe a method allowing establishing of a 
connection between geologic data obtained by coring or logging in wells and seismic data 
obtained by means of a seismic exploration survey, which is based on a statistical 
calibration technique with statistical calibration obtained by bringing together the local 
geologic data (measured in wells) and seismic attributes read on seismic traces obtained in 
the immediate vicinity of each well. 

[0007] The seismic attributes that are subjected to such an interpretation are 

calculated in the post-stack or in the pre-stack domain. 

[0008] In the post-stack domain, the conventional attributes are calculated from the 

amplitudes at the reservoir level, or from the impedance P estimated by means of a 
stratigraphic inversion. 

[0009] In the pre-stack domain, the number of pertinent attributes can increase 

considerably. In this case, the attributes are either the amplitudes at different offsets or 
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incidence angles, or parameters obtained from a pre-stack joint stratigraphic inversion 
(impedances P and S, density, product of the Lame parameters with the density, etc.). 

[0010] Concerning the interpretation of logs, the situation is similar to a very large 

number of measurements being available, hence a large number of attributes allowing 
characterization of the logs, these attributes being often organized in families of equal 
physical significance (for example family of the resistivity attributes, family of the 
attributes associated with the natural or induced radioactivity of the formations, etc.). 
[001 1] A method referred to as principal components analysis (PC A) method, well 

known in the art, can be applied to these attributes in order to analyze the relations existing 
between them and to reduce the number of significant attributes. The principal 
components extracted from the PCA are the new attributes : they define an orthogonal or 
non-orthogonal (rotation) vector base and are linear combinations of initial variables. They 
can be used in statistical pattern recognition algorithms and correspond to a multivariate 
filtering of the initial seismic or log information. Various applications of the PCA method 
are for example described in : 

- Dumay, J., Fournier, F., 1988, "Multivariate Statistical Analyses Applied to Seismic 
Facies Recognition", Geophysics, 53, 1151-1159, 

- Hagen, D., C, 1982, "The Application of Principal Component Analysis to Seismic Data 
Sets", Geoexpl., 20, 93-111, or 

-French Patent A-2,772,138 and corresponding US Patent 6,345,108 filed by the assignee. 
[0012] However, PCA does not take into the attributes group and, consequently, it 

is often difficult to give a physical interpretation of the principal components, or to relate 
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them clearly to the initial attributes, especially if they are organized in groups of equal 
physical significance. 

[0013] The analysis referred to as canonical analysis is also a statistical method 

well known in the art, which allows establishing the relations that may exist between two 
sets of variables in order to know if these two sets describe the same properties. This 
method is for example described by : 

- Hotelling, H., 1936, "Relations Between Two Sets of Variables", Biometrika, 28, 321- 
377. 

[0014] An example of application of canonical analysis is for example described 

by: 

- Founder, F., and Derain, J.F., 1995, "A Statistical Methodology for Deriving Reservoir 
Properties from Seismic Data", Geophysics, 60, 1437-1450. 

[0015] This method is limited to the study of two sets of variables and does 

therefore generally not apply to the multicube seismic information or to the multidomain 
log information. Furthermore, it defines synthetic variables in each subspace associated 
with the two sets, and not a single vector base allowing packing the whole of the initial 
variables and to describe the relation between the two sets. 

[0016] There are several methods allowing generalization of the analysis referred 

to as canonical analysis by extending the analysis to more than two sets of variables. 
Different aspects of the prior art in the field considered are described, for example, in the 
following publications : 
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- Horst, P., 1961, "Relations Among M Sets of Measures", Psychometrika, 26, No.2, 129- 
149, 

- Caroll, J.D., 1968, "A Generalization of Canonical Correlation Analysis to Three or 
More Sets of Variables", Proc. 76 th Conv. Amer. Psych. Ass., 

- Kettenring, J.R., 1971, "Canonical Analysis of Several Sets of Variables", 58, 3, 433- 
450, or 

- Saporta, G., 1990, "Probability, Analyse des Donnees et Statistiques", Technip, Paris. 
[001 7] These methods however have certain limits. In particular, they do not allow 
describing each set separately, the synthetic variables cannot be readily related to the 
various sets and, consequently, be physically interpreted. Furthermore, only the global 
inertia part they represent can be known, and not the inertia part they represent for each 
set. It is therefore not possible to filter each set independently. 

[0018] Another known analysis technique referred to as generalized principal 

components analysis (GPCA) allows comparison of different sets of variables while 
describing each one. It is implemented on data of economic nature for example in : 
- Casin, Ph., 2001, "A Generalization of Principal Component Analysis to K Sets of 
Variables", Computational Statistics & Data Analysis, 35, 417-428. 

SUMMARY OF THE INVENTION 
[0019] The method according to the invention provides an approach wherein, by 

implementing a multivariable statistical technique applied to several sets of seismic 
records or logs available, obtained by exploration of an underground zone, the geophysical 
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events readable on these sets of measurements are packed and filtered. The method 
comprises : 

- distributing these geophysical events into a certain number p of families of geophysical 
attributes having each a particular physical significance; 

- forming, by combination of the geophysical attributes, synthetic variables in much fewer 
number than the number of the attributes, which are obtained by construction of an 
orthogonal vector base in each one of the p analysis sets consisting of the data of each one 
of the families, wherefrom is derived the formation of an orthonormal vector base 
describing the p analysis sets; 

- using the orthonormal vector base for filtering and describing the geophysical events. 
[0020] The methodology for packing and filtering geophysical data in order to 

retrieve therefrom information on the nature of the subsoil is applied to packing and 
filtering sets of multi-domain geophysical data, i.e. either pre-stack seismic data from a 3D 
seismic survey grouped into a certain number of data cubes corresponding each, for 
example, to the same offset (seismic trace cubes referred to as iso-offset cubes) or to the 
same angle of incidence (referred to as iso-incidence angle cubes), or seismic attribute 
cubes obtained by seismic amplitude processing such as, for example, elastic properties 
eubes from a joint stratigraphic inversion of pre-stack seismic data, these cubes being 
obtained by implementing the method described for example in French patent 2,800,473 
(EP-1, 096,270) filed by the assignee. It then shapes the pre-stack seismic information 
contained in seismic trace portions formed by exploration of an underground zone, by 
means of new representative attributes which are thereafter used for geologic 
interpretation or interpretation of an underground reservoir. It can also apply to packing 
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and filtering of pre-stack or post-stack 3D seismic data from repeated surveys (4D seismic 
method) in order to interpret the fluid content variations of a petroleum reservoir under 
production. This method is also applicable to packing of log data of several physical 
natures by forming, from initial measurements, new log attributes in order to improve 
interpretation of the log data in terms of lithologic and petrophysical characteristics of the 
reservoir. 

[0021] Formation of the synthetic variables is obtained for example by combining 

a principal components analysis technique or its derived techniques based on 
diagonalization of a norm (seeking eigenvectors ) and a simple regression technique. 
[0022] According to an implementation mode, the synthetic variables can be 

subjected to a rotation so as to improve the correlations between initial variables, that is 
the geophysical attributes, and synthetic variables. 

[0023] The synthetic variables can be formed from families comprising for 

example one or more seismic trace cubes obtained after a 3D seismic survey, each one 
corresponding to the same offset. 

[0024] Synthetic variables can also be formed from families consisting of seismic 

trace cubes obtained after a 3D seismic survey, each one corresponding to the same angle 
of incidence. 

[0025] Synthetic variables can also be formed from families consisting of seismic 

attribute cubes formed by processing the seismic trace amplitudes obtained at the end of a 
3D seismic survey. 

[0026] Synthetic variables can also be formed from families consisting of seismic 

trace cubes obtained after successive seismic surveys in the zone. 



[0027] It is also possible to form synthetic variables from families consisting of 

seismic attribute cubes formed by processing the seismic trace amplitudes obtained after 
successive seismic surveys in the zone. 

[0028] Synthetic variables can also be formed from families consisting of 

measurements of log attributes of different physical nature, such as the family of 
resistivity measurements, the family of measurements related to the formation 
radioactivity, the family of measurements related to the formation porosity, etc. 
[0029] The method is particularly advantageous in that it allows, when applied to 

seismic or log data : 

- to define new geophysical attributes, linear combinations of the initial geophysical 
attributes, summarizing the initial seismic or log data and allowing filtering of this data, 
while respecting the multidomain structure of the data; 

- to show if all the families of initial geophysical attributes, or only some of them, provide 
redundant information, and to measure the part of this redundant information; and 

- to give a physical interpretation to the new geophysical attributes. 

[0030] The new geophysical attributes thus formed can then be analysed in order to 

retrieve the geologic properties of the seismic or log data. Two classes of methods can be 
conventionally used therefore : a first one using statistical estimation techniques and 
providing quantitative interpretation in terms of geologic properties, and a second one 
using pattern recognition techniques to analyze the seismic facies or the log facies 
(electrofacies) and provide qualitative information. 
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[0031] Examples of these two types of methods are respectively implemented in 

the aforementioned EP Patent 671,017 and US Patent 5,638,269 and French Patent 
2,768,818 and corresponding US Patent 6,051,651 filed by the applicant. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0032] - Figs. 1A to 1C show the organization of the seismic records in p sets : 

each set groups in a cube all of the trace portions corresponding to an identical physical 
criterion; 

[0033] - Figs. 2A to 2C respectively show the maps of the first three new seismic 

attributes Z (1) ,Z (2) and Z (3) ; 

[0034] - Figs. 3 A to 3D respectively show the first true temporal plane (3A) of the 

analysis window and the same plane reconstituted/filtered with 1 (3B), 3 (3C) or 8 (3D) of 
the new seismic attributes for the cube of impedance P; 

[0035] - Figs. 4A to 4D respectively show the first true temporal plane (4A) of the 

analysis window and the same plane reconstituted/filtered with 1 (4B), 3 (4C) or 8 (4D) of 
the new seismic attributes for the cube of impedance S; 

[0036] - Figs. 5A to 5D respectively show the first true temporal plane (5A) and 

the same plane reconstituted/filtered with 1 (5B), 3 (5C) or 8 (5D) of the new seismic 
attributes for the cube of density p; 

[0037] - Figs. 6 respectively show the vertical sections, passing through well W2, 

true (6A), reconstituted/filtered with the first 8 new attributes (6B), and the differences 
(residues) (6C) between these two sections for the cube of impedance P; 
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[0038] - Figs. 7A-C respectively show the vertical sections, passing through well 

W2, true (7A), reconstituted/filtered with the first 8 new attributes (7B), and the 
differences (residues) (7C) between these two sections for the cube of impedance S; 

[0039] - Fig. 8 shows the location of the learning traces used to calibrate the 

classification function within the context of the supervised pattern recognition of the first 8 
new seismic attributes; 

[0040] - Figs. 9A and 9B respectively show the assignments of the seismic events 

and the most reliable assignments of these events within the context of the supervised 
pattern recognition of the first 8 new seismic attributes; and 

[0041] - Figs. 10A and 10B respectively show the assignments of the seismic 

events and the assignments of these seismic events within the context of the non- 
supervised pattern recognition of the first 8 new seismic attributes. 

DETAILED DESCRIPTION OF THE INVENTION 

[0042] Several sets of geophysical records obtained by seismic exploration or 

logging in an underground zone studied are available. The geophysical events to be 
packed/filtered are characterized by geophysical attributes often organized in sets or 
families of particular physical significance (Figures 1A to 1C). These attributes can be the 
succession of the amplitude values of seismic trace portions or of elastic parameters from 
a pre-stack joint elastic inversion, obtained for example by means of the method described 
in the aforementioned French Patent 2,800,473 and EP 1,096,270 along the time window 
studied. These attributes can be, for a given depth, the measurements of a single set of logs 
obtained along a well, such as for example the values of the different resistivity tools, the 
values of the different tools measuring the porosity, etc. 
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1 - Notations 

[0043] Let X be a matrix consisting of p juxtaposed tables Xi , . . . Xj . . . X p : 

X={Xi ...Xj...X p } 9 Xj = [Xj,i ...Xj, k ...Xj, mj ] 

where Xj is a table consisting of mj variables X jsk of n individuals (in the statistical sense 
of the term). Matrix X groups together all the geophysical events available. These events 
are first organized in p sets or families Xj (Figure 1A). Each family Xj corresponds to a 
particular physical property. For analysis of the 3D seismic survey, this family has the 
form of a cube of values corresponding to this property, located by their spatial position 
(line, cdp, time/depth) (Figures IB and 1C). These cubes can for example consist of a set 
of amplitude values of n traces or trace portions associated with a common offset or a 
common angle or the spatial distribution of the values of a common attribute derived from 
these amplitudes such as, for example, the acoustic impedance in P or S waves, the density 
or Lame parameters, in a selected time window of size mj. In this example, each variable 
X jjk groups together the values of the seismic attribute considered (seismic amplitude or 
elastic property) for a temporal sample within the analysis window. In this example, a 
cube comprises as many individuals as traces available. For analysis of the logs, this 
family has the form of a sequence of mj series of values, located by their position along the 
well (time or depth). This sequence can consist, for example, of all the logs of the tools of 
same physical nature, such as the various tools measuring the formation porosity, or the 
tools measuring various resistivity types, on a set of n time or depth points along one or 
more wells at the level of the underground formation studied. In this example, each 
variable X j)k groups together the values of the log attribute considered (resistivity, 
porosity, radioactivity,...) for a time or depth point within the analysis window. In this 
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example, a sequence comprises as many individuals as time or depth points available 
along the well(s) analysed. 

[0044] The m, variables Xj, k of a table Xj define a space denoted by Sj and are 

assumed to be centred. Let D n be the diagonal matrix of the weights of the individuals. 

[0045] Let ITj be an orthogonal projection operator on space Sj defined by : 

Il J =X J (X> J D n X j yX> J D n 

2 - Calculation of the new attributes 

[0046] The purpose of Generalized Principal Components Analysis (GPCA) is 

both to describe the proximities between the p families of variables and to describe 
separately each set by means of synthetic variables that are linear combinations of the 
initial variables. Then, what is common to all the sets or families, or to some of them only, 
and consequently what is different can be shown. GPCA allows reduction of the number 
of variables used to describe a data set and it allows multivariate filtering while keeping 
the multidomain structure of the initial data. To reach this objective, the first synthetic 
variable calculated, called Z (1) , has to be such that: 

1. The p orthogonal projections Zj (1) of Z (l) on spaces Sj have to be as close as possible 
to one another; 

2. In each group j, the dispersion of the n individuals projected on Zj (1) is maximum. 
[0047] This is mathematically translated by the following conditions : 

2 7=1 

ntj _ 

2.m^£Cov 2 {Z?\X hv ) 

z > v=l 
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where R 2 is the multiple correlation coefficient and Zj (1) is the projection of Z ( on space 
Sj normalized at L 

Z (i) =n o) Z o) 

[0048] A new synthetic variable Z (2) is then sought according to the same criteria. 

An orthogonality condition is added between the projections of variables Z (k) of the same 
set. The procedure stops when an orthogonal base is found in each space Sj. The number 
of variables Z (k) does not exceed the largest dimension of spaces Sj. 

[0049] Rotation methods well known to the man skilled in the art can be 

implemented in each space thus found in order to optimize the correlation coefficients 
between the new variables and the initial variables. 

[0050] The first variable Z (1) is the normalized eigenvector of covariance matrix 

XX l D„ corresponding to the greatest eigenvalue. Z (1) is then projected on the p spaces Sj 
and the regression of the my variables by Zj (1) is calculated in each space. The matrix of the 
p juxtaposed residues Xj (2) defines the new space studied. This new space corresponds to 
the initial space from which the influence of Zj (I) has been removed in each set Sj. The mj 
variables of table Xj (2) define a subspace Sj (2) of Sj. The second component Z (2) is the 
normalized eigenvector of covariance matrix X (2) X (2)t D n of the residues corresponding to 
the greatest eigenvalue. Z (2) is projected on the p spaces Sj (2) . The regression of the mj 
variables by Zj (2) is calculated in each space Sj (2) . The regression carried out in each space 
Sj ensures that the projection of Z (2) on space Sj (2) is orthogonal to the projection of Z (1) on 
space Sj. 

[0051] Calculation of the eigenvectors is carried out by means of any method 

known in the art. 
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[0052] In practice, matrix XX l D n being of order n, where n can be very great, the 

normalized eigenvector V (k) of matrix X^DnX 00 associated with the greatest eigenvalue is 
calculated in a stage k. Vector Z (k) is then simply calculated by linear transformation Z 00 = 
x (k) v (k) and normalized. 

[0053] Each set can be reconstituted with a small number of attributes. The initial 

variables of set j can be expressed as the sum of a linear combination of variables Zj (k) and 
of a residue. If all the variables Zj (k) of a set are used, the residue is zero. The coefficients 
of the linear combination are the coefficients calculated after each regression in spaces Sj. 

[0054] One of the properties of the method according to the invention is that 

variables Z (k) are orthogonal two by two. 

[0055] For interpretation of the synthetic variables calculated, then calculate for 

each set j: 

1. The square of the correlation between variable Z (k) and variable Zj (k) of set j, 
indicating the degree of proximity between these two variables; and 

2. the inertia percentage of set j explained by variable Zj (k) . 

[0056] The correlations between variables Z (k) and initial variables X jjk are also 

calculated to show possible relations between the X jjk . The correlations between variables 
Zj (k) and Xj, k of the same set j, or the correlations between variables Zj (k) for the same k 
can also be calculated. 

[0057] The new synthetic variables thus created constitute the new geophysical 

attributes sought. 

3 - Application example 
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[0058] The method according to the invention was applied to three cubes of elastic 

properties (impedance P denoted by Ip, impedance S denoted by Is and density denoted by 
p) obtained after pre-stack joint elastic inversion, by means of the method described in the 
aforementioned French Patent 2,800,473, of seismic records grouped in cubes of common 
incidence angles. Each elastic property generates a set of seismic attributes or variables. 
The analysis time window consists of 19 signal samples taken from seismic records. The 
values of the elastic property for each temporal sample within the analysis window to be 
variables are considered. Each set thus consists of 19 variables made up of as many 
individuals as there are traces available. Since there are three variable sets, a total 57 
variables exists. The three parameter cubes have been flattened in relation to a reference 
surface. 

[0059] The synthetic variables Z (k) obtained summarize the respective information 

of the three elastic properties, while taking into account the relations between them. 

3.1 - Analysis of the relations between impedances P, S and density p 
[0060] Table 1 hereunder gives the variance proportion in percentage explained by 

the projections of the first 8 synthetic variables Z (1) to Z (8) in the three sets. It shows that 8 
variables are sufficient to summarize 87.4 % of the variance of impedance P, 94.4 % of the 
variance of impedance S and 88.4 % of the variance of the density. 
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[0061] Table 2 hereafter gives the square of the correlation between synthetic 

variables Z 00 and their projections on the three sets. This value indicates the degree of 
proximity between the variable and the set. The closer the value is to 1, the closer the 
variable is to its projection and therefore to the set. 
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[0062] For example, Z (1) is very close to Is and close to p, but far from Ip, whereas 

Z (2) is close to Ip but very far from Is and p. Out of the first 8 variables Z (k) , three are close 
to Ip alone, two are close to both Is and p, one is close to Is alone, and two are close to p 
alone. It thus appears that Ip and Is do not contain the same information and that part of 
the information contained in Is and in p seems to be common. 

3.2 - Filtering of impedances P, S and of the density 

[0063] The three elastic properties can be reconstituted with a limited number of 

new seismic attributes. The number of synthetic variables used for reconstitution can be 
different from one set to the next : it depends on the inertia percentage to be 
reconstituted/filtered. 

3.3 - Supervised and non-supervised pattern recognition of the new attributes 
[0064] It is possible to interpret geologically the new seismic attributes through a 
facies analysis with a supervised or non-supervised pattern recognition algorithm. 
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[0065] In supervised pattern recognition, it is necessary to select learning classes 

grouping together seismic events located in the vicinity of wells whose geologic 
interpretation has been carried out. This methodology allows classification of all of the 
seismic events according to the geologic characters observed in the wells. 

[0066] In the present case, the geologic knowledge of the three available wells Wl, 

W2 and W3, and the general geologic knowledge of the zone studied lead us to 
considering three classes : a first class corresponding to wells Wl and W2 represents 
turbidite sands, a second class corresponding to well W3 represents turbidite muds and a 
third class represents hemipelagites. Wells Wl and W2 showing close geologic 
characteristics, only the seismic events located around well W2 are used, well Wl serving 
as a blind test. The location of well W3 being uncertain, the seismic events used to 
characterize class 2 are located around a virtual well defined in the proximity of W3. 
Finally, since there are no well for class 3, a virtual well is defined from our general 
geologic knowledge. Each class consists of 121 seismic events (llxl 1). These events are 
used to calibrate a classification function which is then applied to the whole of the seismic 
events. 

[0067] Figs. 8A and 8B respectively show the assignments obtained and the most 

probable assignments, i.e. having an assignment probability greater than or equal to 80 %. 
Well Wl, used in blind test, is correctly assigned to the same class as well W2. Well W3 is 
located at the boundary of classes 1 and 2. Class 2 is distributed along a channeling form. 
The map of Fig. 8B shows that the results are stable : few events exhibit a bad assignment 
probability. 
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[0068] In non-supervised pattern recognition, the learning classes are formed using 

the modes of a multivariate probability density function calculated from the whole of the 
seismic events, in a space generated by the first eight new attributes. The seismic events of 
higher probability density that constitute the peaks of this density function are selected to 
carry out learning of the classification function. This function is then applied to the whole 
of the seismic events considered to assign them to one of the classes : figs. 9A and 9B 
show the assignments obtained, as well as the most probable assignments, i.e. having an 
assignment probability greater than or equal to 80 %. The maps obtained are very close to 
the maps of Figs. 8A and 8B. Class 3 of the supervised analysis corresponds to classes 3, 4 
and 5 of the non-supervised analysis. Wells Wl and W2 are assigned to the same class and 
well W3 is still located at the boundary between classes 3, and especially 1 and 2. This 
convergence of results validates the interpretation of the supervised approach. 
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